home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PMAT12 / PMAT.PAS < prev    next >
Pascal/Delphi Source File  |  1993-01-23  |  24KB  |  895 lines
  1. {
  2. P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
  3.  
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/30/93
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9.  
  10.  
  11. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  12. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  13. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  14. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  15. FROM USE OF THIS PROGRAM.
  16. }
  17.  
  18. Unit pmat;
  19.  
  20. Interface
  21.  
  22. Type
  23.      fp = ^double;
  24.      ap = Array[1..1] Of double;
  25.      apptr = ^ap;
  26.      app = Array[1..1] Of apptr;
  27.  
  28.      vmatrixptr = ^vmatrix;
  29.      vmatrix = Object
  30.         r, c : integer;
  31.         Function m( i, j: integer ): double;
  32.         Function mm( i, j: integer ) : fp;
  33.         constructor MakeMatrix( vr, vc: integer );
  34.         Destructor KillVmatrix;
  35.         Procedure Garbage;
  36.         Procedure Show( strng: String );
  37.         Procedure infomatrix( strng: String );
  38.  
  39.         private
  40.            v,vcheck: ^app;
  41.            nelems : longint;
  42.            onstack : boolean;
  43.            Procedure purgevectors;
  44.            Procedure allocvectors( rr, cc: integer );
  45.      End;
  46.  
  47.      mstackptr = ^mstack;
  48.      mstack = Object
  49.         constructor InitMatrixStack;
  50.         Destructor KillMatrixStack;
  51.         Procedure inclevel;
  52.         Procedure declevel;
  53.         Function ReturnMat: vmatrixptr;
  54.         Function DecReturn: vmatrixptr;
  55.         Procedure push( Var a: vmatrixptr );
  56.         Procedure nrerror( strng: String );
  57.         Procedure DumpStack;
  58.  
  59.         private
  60.           nummats, level : integer;
  61.           next : mstackptr;
  62.           mv : vmatrixptr;
  63.           Procedure cleanstack( a: vmatrixptr );
  64.           Function pop: vmatrixptr;
  65.      End;
  66.  
  67. Var
  68.     dispatch : mstackptr;
  69.  
  70. Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) : vmatrixptr;
  71. Function ident( n: integer ): vmatrixptr;
  72.  
  73. Function Mult( A, B: vmatrixptr ):vmatrixptr;
  74. Function emult( a, b: vmatrixptr ):vmatrixptr;
  75. Function scmult( a: double; b: vmatrixptr ): vmatrixptr;
  76. Function multsc( b: vmatrixptr; a: double ): vmatrixptr;
  77.  
  78. Function divis( a, b: vmatrixptr ):vmatrixptr;
  79. Function scdivis( a: double; b: vmatrixptr ): vmatrixptr;
  80. Function divissc( b: vmatrixptr; a: double ): vmatrixptr;
  81.  
  82. Function add( a, b: vmatrixptr ): vmatrixptr;
  83. Function scadd( a: double; b: vmatrixptr ): vmatrixptr;
  84. Function addsc( b: vmatrixptr; a: double ): vmatrixptr;
  85.  
  86. Function neg( A: vmatrixptr ): vmatrixptr;
  87. Function sub( A, B: vmatrixptr ):vmatrixptr;
  88. Function scsub( A: double; B: vmatrixptr ):vmatrixptr;
  89. Function subsc( B: vmatrixptr; A: double ):vmatrixptr;
  90.  
  91. Function tran( a: vmatrixptr ): vmatrixptr;
  92. Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr;
  93. Function inv( Amat: vmatrixptr ): vmatrixptr;
  94.  
  95. Function fill( rr, cc: integer; d: double ):vmatrixptr;
  96. Function submat( a: vmatrixptr; lr, r, lc, c: integer ): vmatrixptr;
  97. Function cv( a, b: vmatrixptr ): vmatrixptr;
  98. Function ch( a, b: vmatrixptr ): vmatrixptr;
  99. Function vecdiag( a: vmatrixptr ): vmatrixptr;
  100.  
  101. Function reada( fid: String ): vmatrixptr;
  102. Procedure writea( fid: String; a: vmatrixptr; titlestr: String );
  103.  
  104. Procedure setwid( wid: integer );
  105. Procedure setdec( decimals: integer );
  106. Function getwid: integer;
  107. Function getdec: integer;
  108.  
  109. Implementation
  110.  
  111. {///////////////// vmatrix functions ///////////////}
  112. {$R-}
  113. Procedure vmatrix.allocvectors;
  114. Var i: integer;
  115. Begin
  116.     if c > 8192 then 
  117.        dispatch^.nrerror('A matrix cannot have more than 8192 columns');
  118.     getmem( v, r * sizeof( apptr ) );
  119.     For i := 1 To r Do
  120.         getmem( v^[i], c * sizeof( double ) );
  121. End;
  122.  
  123. Procedure vmatrix.purgevectors;
  124. Var i: integer;
  125. Begin
  126.     If v = Nil Then exit;
  127.     For i := r Downto 1 Do
  128.         freemem( v^[i], c * sizeof( double ) );
  129.     freemem( v, r * sizeof( apptr ) );
  130. End;
  131. {$R+}
  132.  
  133. constructor vmatrix.MakeMatrix;
  134. Begin
  135.     r := vr;
  136.     c := vc;
  137.     onstack := false;
  138.     if c > 8192 then 
  139.        dispatch^.nrerror('A matrix cannot have more than 8192 columns');
  140.     nelems := longint( longint( r ) * longint( c ) ) * longint( sizeof( double ) );
  141.     If nelems < maxAvail - 256 Then allocvectors( r, c )
  142.     Else dispatch^.nrerror( 'cannot allocate matrix' );
  143.     vcheck := v;
  144. End; { MakeMatrix }
  145.  
  146. Destructor vmatrix.KillVmatrix;
  147. Begin
  148.     purgevectors;
  149.     vcheck := Nil;
  150. End;
  151.  
  152. Procedure vmatrix.Garbage;
  153. Var errorint : boolean;
  154. Begin
  155.     errorint := false;
  156.     If vcheck <> v Then errorint := true;
  157.     If v = Nil Then errorint := true;
  158.     If r < 1 Then errorint := true;
  159.     If c < 1 Then errorint := true;
  160.     If errorint Then dispatch^.Nrerror( 'garbage' );
  161. End;
  162.  
  163. {$R-}
  164.  
  165. Function vmatrix.m;
  166. Begin
  167.     { access a member on the right side of assignment }
  168.     If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
  169.         dispatch^.nrerror( 'm: index out of range' );
  170.     m := v^[i]^[j];
  171. End;
  172.  
  173. Function vmatrix.mm;
  174. Begin
  175.     { store a matrix element on the left side of assignment }
  176.     If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
  177.         dispatch^.nrerror( 'mm: index out of range' );
  178.     mm := @v^[i]^[j];
  179. End;
  180.  
  181. {$R+}
  182.  
  183. Procedure CopySD( Var source, dest: vmatrixptr );
  184. Var
  185.   i,j : integer;
  186. Begin
  187.     source^.Garbage;
  188.     dest^.Garbage;
  189.     
  190.     If source = dest Then exit;
  191.     If source^.onstack Then
  192.         With dest^ Do Begin 
  193.             purgevectors;
  194.             nelems := source^.nelems;
  195.             v := dispatch^.Pop^.v;
  196.             r := source^.r;
  197.             c := source^.c;
  198.             vcheck := v;
  199.             source^.v := Nil;
  200.             dispose( source, killvmatrix );
  201.             exit;
  202.         End;
  203.     
  204.     If source^.v = dest^.v Then
  205.         With dest^ Do Begin 
  206.             r := source^.r;
  207.             c := source^.c;
  208.             allocvectors( r, c );
  209.             vcheck := v;
  210.             nelems := source^.nelems;
  211.         End;
  212.     
  213.     If ( dest^.r <> source^.r ) Or( dest^.c <> source^.c ) Then
  214.         With dest^ Do Begin 
  215.             purgevectors;
  216.             r := source^.r;
  217.             c := source^.c;
  218.             allocvectors( r, c );
  219.             vcheck := v;
  220.             nelems := source^.nelems;
  221.         End;
  222.     
  223.     For i := 1 To dest^.r Do
  224.         For j := 1 To dest^.c Do
  225.             dest^.mm( i, j )^ := source^.m( i, j );
  226. End;
  227.  
  228. Function matequals;
  229.      Begin
  230.          rop^.Garbage;
  231.          lop^.Garbage;
  232.          copySD( rop, lop );
  233.          dispatch^.CleanStack( lop );
  234.          matequals := lop;
  235.      End;
  236.  
  237.  
  238. {////////////////////// stack functions /////////}
  239.  
  240.  
  241. constructor mstack.InitMatrixStack;
  242. Begin
  243.     nummats := 0;
  244.     level := 1;
  245.     mv := Nil;
  246.     getmem( next, sizeof( mstack ) );
  247.     With next^ Do Begin 
  248.         level := 0;
  249.         next := Nil;
  250.         mv := Nil;
  251.         nummats := 0;
  252.     End;
  253. End;
  254.  
  255. Destructor mstack.KillMatrixStack;
  256. Begin
  257.     freemem( next, sizeof( mstack ) );
  258. End;
  259.  
  260. Procedure mstack.nrerror;
  261. Begin
  262.     writeln( 'Fatal Error: ', strng );
  263.     writeln( 'Exiting to system' );
  264.     halt;
  265. End;
  266.  
  267. Procedure mstack.inclevel;
  268. Begin
  269.     level := level + 1;
  270. End;
  271.  
  272. Procedure mstack.declevel;
  273. Begin
  274.     level := level - 1;
  275.     If level <= 0 Then nrerror( 'declevel: decremented too often' );
  276. End;
  277.  
  278. Function mstack.ReturnMat;
  279. Begin
  280.     ReturnMat := next^.mv;
  281. End;
  282.  
  283. Function mstack.DecReturn;
  284. Begin
  285.     declevel;
  286.     DecReturn := next^.mv;
  287. End;
  288.  
  289. Procedure mstack.push;
  290.    Var
  291.       t : mstackptr;
  292. Begin
  293.     a^.Garbage;
  294.     getmem( t, sizeof( mstack ) );
  295.     t^.level := dispatch^.level;
  296.     t^.next := dispatch^.next;
  297.     t^.mv := a;
  298.     a^.onstack := true;
  299.     dispatch^.next := t;
  300.     dispatch^.nummats := dispatch^.nummats + 1;
  301.     t^.nummats := dispatch^.nummats;
  302. End;
  303.  
  304. Function mstack.pop;
  305. Var t : mstackptr;
  306.     vv : vmatrixptr;
  307. Begin
  308.     If dispatch^.level < 0 Then
  309.         nrerror( 'pop: missmatched inclevel-declevel statments, level < 0' );
  310.     If ( dispatch^.next^.next = Nil ) Or( Nil = dispatch^.next^.mv ) Then
  311.         nrerror( 'pop: popping off the tail' );
  312.     
  313.     t := dispatch^.next;
  314.     dispatch^.next := dispatch^.next^.next;
  315.     dispatch^.nummats := dispatch^.nummats - 1;
  316.     t^.mv^.Garbage;
  317.     vv := t^.mv;
  318.     freemem( t, sizeof( mstack ) );
  319.     t := Nil;
  320.     vv^.onstack := false;
  321.     pop := vv;
  322. End;
  323.  
  324. Procedure mstack.cleanstack;
  325. Var vv : vmatrixptr;
  326. Begin
  327.     While ( dispatch^.next^.next <> Nil ) And
  328.             ( dispatch^.next^.level >= dispatch^.level ) Do Begin 
  329.         vv := pop;
  330.         If vv <> a Then
  331.             dispose( vv, KillVmatrix );
  332.     End;
  333. End;
  334.  
  335. Procedure mstack.dumpstack;
  336. Var mm : mstackptr;
  337. Begin
  338.     mm := dispatch;
  339.     writeln( '------------ Stack Dump --------------' );
  340.     While mm <> Nil Do Begin 
  341.         With mm^ Do Begin 
  342.             writeln( '----------------------------------------------' );
  343.             writeln( 'nummats, level, next, vector: ', nummats, ' ', level, ' ', 
  344.                 longint( next ), ' ', longint( mv ) );
  345.             If mv <> Nil Then mv^.infomatrix( 'stack matrix information' )
  346.             Else
  347.                 writeln;
  348.         End;
  349.         mm := mm^.next;
  350.     End;
  351. End;
  352. {////////////////// matrix functions ///////////////////}
  353.  
  354. Procedure vmatrix.Show;
  355.    Var i,j,w,d: integer;
  356.    Begin
  357.        w := getwid;
  358.        d := getdec;
  359.        writeln;
  360.        writeln( strng );
  361.        writeln;
  362.        For i := 1 To r Do Begin 
  363.            For j := 1 To c Do
  364.                write( m( i, j ): w: d, ' ' );
  365.            writeln;
  366.        End;
  367.        writeln;
  368.    End; { show }
  369.  
  370. Procedure vmatrix.infomatrix;
  371.    Begin
  372.        writeln;
  373.        writeln( strng );
  374.        writeln;
  375.        writeln( 'r,c,size : ', r, ' ', c, ' ', nelems, ' bytes' );
  376.        writeln( 'v, vcheck: ', longint( v ), ' ', longint( vcheck ) );
  377.    End; { show }
  378.  
  379. Function mult;
  380.     Var i,j,k: integer;
  381.         C: vmatrixptr;
  382.         s: extended;
  383.     Begin
  384.         a^.Garbage;
  385.         b^.Garbage;
  386.         If A^.c <> B^.r Then
  387.             dispatch^.nrerror( 'mult: matrices do not conform' )
  388.         Else Begin 
  389.             new( C, MakeMatrix( a^.r, b^.c ) );
  390.             For i := 1 To A^.r Do
  391.                 For j := 1 To B^.c Do Begin 
  392.                     s := 0;
  393.                     For k := 1 To B^.r Do
  394.                         s := s + A^.m( i, k ) * B^.m( k, j );
  395.                     C^.mm( i, j )^ := s;
  396.                 End;
  397.         End;
  398.         Dispatch^.Push( C );
  399.         Mult := Dispatch^.ReturnMat;
  400.     End; {mult}
  401.  
  402. Function Emult{( A, B: vmatrixptr): vmatrixptr};
  403.     Var i,j: integer;
  404.         C: vmatrixptr;
  405.     Begin
  406.         a^.Garbage;
  407.         b^.Garbage;
  408.         If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
  409.             dispatch^.Nrerror( 'Emult: matrices do not conform' )
  410.         Else Begin 
  411.             new( c, MakeMatrix( A^.r, A^.c ) );
  412.             For i := 1 To a^.r Do
  413.                 For j := 1 To a^.c Do
  414.                     C^.mm( i, j )^ := A^.m( i, j ) * B^.m( i, j );
  415.         End;
  416.         Dispatch^.push( c );
  417.         emult := Dispatch^.returnmat;
  418.     End; {emult}
  419.  
  420. Function scmult{( a: double; b: vmatrixptr): vmatrixptr};
  421.     Var i,j: integer;
  422.         C: vmatrixptr;
  423.     Begin
  424.         b^.Garbage;
  425.         new( c, MakeMatrix( b^.r, b^.c ) );
  426.         For i := 1 To b^.r Do
  427.             For j := 1 To b^.c Do
  428.                 C^.mm( i, j )^ := a * B^.m( i, j );
  429.         Dispatch^.push( c );
  430.         scmult := Dispatch^.returnmat;
  431.     End; {scmult}
  432.  
  433. Function multsc{( b: vmatrixptr; a: double): vmatrixptr};
  434.     Var i,j: integer;
  435.         C: vmatrixptr;
  436.     Begin
  437.         b^.Garbage;
  438.         new( c, MakeMatrix( b^.r, b^.c ) );
  439.         For i := 1 To b^.r Do
  440.             For j := 1 To b^.c Do
  441.                 C^.mm( i, j )^ := a * B^.m( i, j );
  442.         Dispatch^.push( c );
  443.         multsc := Dispatch^.returnmat;
  444.     End; {multsc}
  445.  
  446. Function divis{( A, B: vmatrixptr): vmatrixptr};
  447.     Var i,j: integer;
  448.         C: vmatrixptr;
  449.         dd: double;
  450.     Begin
  451.         a^.Garbage;
  452.         b^.Garbage;
  453.         If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
  454.             dispatch^.Nrerror( 'Emult: matrices do not conform' )
  455.         Else Begin 
  456.             new( c, MakeMatrix( b^.r, b^.c ) );
  457.             For i := 1 To a^.r Do
  458.                 For j := 1 To a^.c Do Begin 
  459.                     dd := b^.m( i, j );
  460.                     If dd = 0.0 Then
  461.                         dispatch^.nrerror( 'DIVIS: zero divisor' );
  462.                     C^.mm( i, j )^ := A^.m( i, j ) / dd;
  463.                 End;
  464.         End;
  465.         Dispatch^.push( c );
  466.         divis := Dispatch^.returnmat;
  467.     End; {divis}
  468.  
  469. Function scdivis{( A : double; B: vmatrixptr): vmatrixptr};
  470.     Var i,j: integer;
  471.         C: vmatrixptr;
  472.         dd: double;
  473.     Begin
  474.         b^.Garbage;
  475.         new( c, MakeMatrix( b^.r, b^.c ) );
  476.         For i := 1 To b^.r Do
  477.             For j := 1 To b^.c Do Begin 
  478.                 dd := b^.m( i, j );
  479.                 If dd = 0.0 Then
  480.                     dispatch^.nrerror( 'SCDIVIS: zero divisor' );
  481.                 C^.mm( i, j )^ := A / dd;
  482.             End;
  483.         Dispatch^.push( c );
  484.         scdivis := Dispatch^.returnmat;
  485.     End; {divis}
  486.  
  487. Function divissc{( A : double; B: vmatrixptr): vmatrixptr};
  488.     Var i,j: integer;
  489.         C: vmatrixptr;
  490.     Begin
  491.         b^.Garbage;
  492.         If a = 0.0 Then
  493.             dispatch^.nrerror( 'DIVISSC: zero divisor' );
  494.         new( c, MakeMatrix( b^.r, b^.c ) );
  495.         For i := 1 To b^.r Do
  496.             For j := 1 To b^.c Do Begin 
  497.                 C^.mm( i, j )^ := b^.m( i, j ) / A;
  498.             End;
  499.         Dispatch^.push( c );
  500.         divissc := Dispatch^.returnmat;
  501.     End; {divissc}
  502.  
  503. Function ADD{( A, B: vmatrixptr): vmatrixptr};
  504.     Var i,j: integer;
  505.         C: vmatrixptr;
  506. Begin
  507.     a^.Garbage;
  508.     b^.Garbage;
  509.     If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
  510.         dispatch^.Nrerror( 'add: matrices do not conform' )
  511.     Else Begin 
  512.         new( c, MakeMatrix( A^.r, A^.c ) );
  513.         For i := 1 To a^.r Do
  514.             For j := 1 To a^.c Do
  515.                 C^.mm( i, j )^ := A^.m( i, j ) + B^.m( i, j );
  516.     End;
  517.     Dispatch^.push( c );
  518.     ADD := Dispatch^.returnmat;
  519. End; {add}
  520.  
  521. Function scADD{( A : double; B: vmatrixptr): vmatrixptr};
  522.     Var i,j: integer;
  523.         C: vmatrixptr;
  524. Begin
  525.     b^.Garbage;
  526.     new( c, MakeMatrix( b^.r, b^.c ) );
  527.     For i := 1 To b^.r Do
  528.         For j := 1 To b^.c Do
  529.             C^.mm( i, j )^ := A + B^.m( i, j );
  530.     Dispatch^.push( c );
  531.     scADD := Dispatch^.returnmat;
  532. End; {scadd}
  533.  
  534. Function ADDSC{(  B : vmatrixptr; A : double): vmatrixptr};
  535.     Var i,j: integer;
  536.         C: vmatrixptr;
  537. Begin
  538.     b^.Garbage;
  539.     new( c, MakeMatrix( b^.r, b^.c ) );
  540.     For i := 1 To b^.r Do
  541.         For j := 1 To b^.c Do
  542.             C^.mm( i, j )^ := A + B^.m( i, j );
  543.     Dispatch^.push( c );
  544.     ADDsc := Dispatch^.returnmat;
  545. End; {addsc}
  546.  
  547. Function neg;
  548.   Var i,j: integer;
  549.   B: vmatrixptr;
  550.   Begin
  551.       a^.garbage;
  552.       new( b, MakeMatrix( a^.r, a^.c ) );
  553.       For i := 1 To a^.r Do
  554.           For j := 1 To a^.c Do
  555.               B^.mm( i, j )^ := - A^.m( i, j );
  556.       dispatch^.Push( b );
  557.       neg := dispatch^.returnmat;
  558.   End; {negate}
  559. Function sub{( A, B: vmatrixptr): vmatrixptr};
  560.     Var i,j: integer;
  561.         C: vmatrixptr;
  562. Begin
  563.     a^.Garbage;
  564.     b^.Garbage;
  565.     If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
  566.         dispatch^.Nrerror( 'sub: matrices do not conform' )
  567.     Else Begin 
  568.         new( c, MakeMatrix( A^.r, A^.c ) );
  569.         For i := 1 To a^.r Do
  570.             For j := 1 To a^.c Do
  571.                 C^.mm( i, j )^ := A^.m( i, j ) - B^.m( i, j );
  572.     End;
  573.     Dispatch^.push( c );
  574.     sub := Dispatch^.returnmat;
  575. End; {sub}
  576.  
  577. Function scsub{ (a: double; b:vmatrixptr): vmatrixptr;};
  578.     Var i,j: integer;
  579.         C: vmatrixptr;
  580.     Begin
  581.         b^.Garbage;
  582.         new( c, MakeMatrix( b^.r, b^.c ) );
  583.         For i := 1 To b^.r Do
  584.             For j := 1 To b^.c Do
  585.                 C^.mm( i, j )^ := A - B^.m( i, j );
  586.         dispatch^.Push( c );
  587.         scSub := dispatch^.returnmat;
  588.     End; { scsub }
  589.  
  590. Function subsc{(b: vmatrixptr; a : double): vmatrixptr;};
  591.     Var i,j: integer;
  592.         C: vmatrixptr;
  593.     Begin
  594.         b^.Garbage;
  595.         new( c, MakeMatrix( b^.r, b^.c ) );
  596.         For i := 1 To b^.r Do
  597.             For j := 1 To b^.c Do
  598.                 C^.mm( i, j )^ := B^.m( i, j ) - A;
  599.         dispatch^.Push( c );
  600.         Subsc := dispatch^.returnmat;
  601.     End; { subsc }
  602.  
  603. Function tran;
  604.     Var i,j: integer;
  605.         c : vmatrixptr;
  606. Begin
  607.     a^.garbage;
  608.     new( c, MakeMatrix( A^.c, A^.r ) );
  609.     For i := 1 To A^.r Do
  610.         For j := 1 To A^.c Do
  611.             c^.mm( j, i )^ := A^.m( i, j );
  612.     dispatch^.Push( c );
  613.     tran := dispatch^.returnmat;
  614. End; {transpose}
  615.  
  616.  
  617. Function sweep;
  618. {
  619.      SUBROUTINE TO SWEEP THE NXN MATRIX A ON ITS K1 THRU K2
  620.      DIACONAL ELEMENTS (SWP(K)SWP(K)A=A)
  621.  
  622.      INPUT : A Matrix to be swept
  623.              K1,K2: elements to sweep; if k1=1 and k2= n then gets inverse
  624. }
  625. Const Eps = 1.0E-8;
  626. Var I,J,k,n,it: integer;
  627.      d,ss: extended;
  628. Begin
  629.     a^.garbage;
  630.     If A^.r <> A^.c Then
  631.         dispatch^.nrerror( 'sweep: Amat not square in Sweep: ' );
  632.     
  633.     n := A^.r;
  634.     If K2 < K1 Then Begin 
  635.         k := k1; k1 := k2; k2 := k;
  636.     End;
  637.     For K := K1 To K2 Do Begin 
  638.         If Abs( A^.m( K, K ) ) < eps Then Begin 
  639.             For it := 1 To n Do Begin 
  640.                 a^.mm( it, k )^ := 0;
  641.                 a^.mm( k, it )^ := 0;
  642.             End;
  643.         End
  644.         Else Begin 
  645.             D := 1.0 / A^.m( K, K );
  646.             A^.mm( K, K )^ := 1.0;
  647.             For I := 1 To N Do A^.mm( K, I )^ := D * A^.m( K, i );
  648.             For J := 1 To N Do
  649.                 If J <> K Then A^.mm( J, k )^ :=
  650.                         - D * A^.m( J, K );
  651.             For J := 1 To N Do
  652.                 If J <> K Then
  653.                     For I := 1 To N Do Begin 
  654.                         ss := A^.m( J, I );
  655.                         If I <> K Then A^.mm( J, I )^ :=
  656.                                 ss + (1.0 / D) * A^.m( J, K ) * A^.m( K, I );
  657.                     End;
  658.         End;
  659.     End;
  660.     dispatch^.push( a );
  661.     sweep := dispatch^.returnmat;
  662. End; {sweep}
  663.  
  664. Function inv;
  665.    Var AMat1: vmatrixptr;
  666.    Begin
  667.        dispatch^.inclevel;
  668.        amat^.garbage;
  669.        If AMat^.r <> AMat^.c Then
  670.            dispatch^.nrerror( 'Matrix not square in Invert ' );
  671.        new( amat1, Makematrix( amat^.r, amat^.r ) );
  672.        amat1 := matequals( amat1, amat );
  673.        amat1^.show( 'amat1' );
  674.        AMat1 := matequals( amat1, Sweep( amat1, 1, amat1^.r ) );
  675.        dispatch^.push( amat1 );
  676.        inv := dispatch^.Decreturn;
  677.    End;  {Inv}
  678.  
  679. {//////////////////// patterned matrices ///////////////////}
  680.  
  681. Function ident;
  682.   Var i,j: integer;
  683.       a: vmatrixptr;
  684. Begin
  685.     new( a, MakeMatrix( n, n ) );
  686.     For i := 1 To n Do
  687.         For j := 1 To n Do
  688.             If i = j Then a^.mm( i, i )^ := 1.0
  689.             Else a^.mm( i, j )^ := 0;
  690.     dispatch^.Push( a );
  691.     ident := dispatch^.ReturnMat;
  692. End; { ident }
  693.  
  694.  
  695. Function fill;
  696.   Var i,j: integer;
  697.       c : vmatrixptr;
  698.   Begin
  699.       If ( rr < 1 ) Or( cc < 1 ) Then
  700.           dispatch^.nrerror( 'fill: negative rows or columns' );
  701.       new( c, makematrix( rr, cc ) )    ;
  702.       For i := 1 To rr Do
  703.           For j := 1 To cc Do
  704.               c^.mm( i, j )^ := d;
  705.       dispatch^.push( c );
  706.       fill := dispatch^.returnmat;
  707.   End;
  708.  
  709. Function submat;
  710.    Var i,j,rc,cc: integer;
  711.        b : vmatrixptr;
  712. Begin
  713.     a^.Garbage;
  714.     If ( lr < 1 ) Or( r > a^.r ) Or( lc < 1 ) Or( c > a^.c ) Then dispatch^.nrerror( 'SUBMAT: index out of range' );
  715.     If ( lr > r ) Or( lc > c ) Then
  716.         dispatch^.nrerror( 'SUBMAT: indexes out of order, lc>c or lr>r' );
  717.     
  718.     new( b, makematrix( r - lr + 1, c - lc + 1 ) );
  719.     
  720.     rc := b^.r;
  721.     cc := b^.c;
  722.     For i := 1 To rc Do
  723.         For j := 1 To cc Do
  724.             b^.mm( i, j )^ := a^.m( lr - 1 + i, lc - 1 + j );
  725.     
  726.     dispatch^.push( b );
  727.     submat := dispatch^.ReturnMat;
  728. End;
  729.  
  730. Function ch;
  731.   Var i,j,k: integer;
  732.       c : vmatrixptr;
  733. Begin
  734.     { horizontal concatination }
  735.     a^.Garbage;
  736.     b^.Garbage;
  737.     If a^.r <> b^.r Then
  738.         dispatch^.nrerror( 'CH: matrices have different number of rows' );
  739.     
  740.     k := a^.c + b^.c;
  741.     new( c, makematrix( a^.r, k ) );
  742.     
  743.     For i := 1 To a^.r Do
  744.         For j := 1 To k Do
  745.             If j <= a^.c Then c^.mm( i, j )^ := a^.m( i, j )
  746.             Else c^.mm( i, j )^ := b^.m( i, j - a^.c );
  747.     
  748.     dispatch^.push( c );
  749.     ch := dispatch^.ReturnMat;
  750. End;     { ch }
  751.  
  752. Function cv;
  753.   Var i,j,k : integer;
  754.           c : vmatrixptr;
  755. Begin
  756.     { vertical concatination }
  757.     a^.Garbage;
  758.     b^.Garbage;
  759.     If a^.c <> b^.c Then
  760.         dispatch^.nrerror( 'CV: matrices have different number of rows' );
  761.     
  762.     k := a^.r + b^.r;
  763.     new( c, MakeMatrix( k, a^.c ) );
  764.     
  765.     For i := 1 To a^.c Do
  766.         For j := 1 To k Do
  767.             If j <= a^.r Then c^.mm( j, i )^ := a^.m( j, i )
  768.             Else c^.mm( j, i )^ := b^.m( j - a^.r, i );
  769.     
  770.     dispatch^.push( c );
  771.     cv := dispatch^.ReturnMat;
  772. End;     { cv }
  773.  
  774. Function vecdiag;
  775. Var  i,j,r: integer;
  776.      b : vmatrixptr;
  777.      t : double;
  778. Begin
  779.     { insert a vector into the diagonal of I }
  780.     a^.Garbage;
  781.     r := 0;
  782.     If a^.r = 1 Then r := a^.c;
  783.     If a^.c = 1 Then r := a^.r;
  784.     If r = 0 Then dispatch^.nrerror( 'vecdiag: a is not a vector' );
  785.     new( b, MakeMatrix( r, r ) );
  786.     
  787.     For i := 1 To r Do Begin 
  788.         If a^.r = 1 Then t := a^.m( 1, i )
  789.         Else  t := a^.m( i, 1 );
  790.         For j := 1 To r Do
  791.             If i = j Then b^.mm( i, j )^ := t
  792.             Else b^.mm( i, j )^ := 0.0
  793.     End;
  794.     
  795.     dispatch^.push( b );
  796.     vecdiag := dispatch^.ReturnMat;
  797. End;     {vecdiag}
  798.  
  799. Function FileExists( FileName: String )
  800.                                 : Boolean;
  801. { Returns True if file exists; otherwise,
  802.   it returns False. Closes the file if
  803.   it exists. }
  804. Var
  805.   f: File;
  806. Begin
  807.     {$I-}
  808.     Assign( f, FileName );
  809.     Reset( f );
  810.     Close( f );
  811.     {$I+}
  812.     FileExists := (IOResult = 0) And
  813.         (FileName <> '');
  814. End;  { FileExists }
  815.  
  816. Function ReadA;
  817. Var i,j  : integer;
  818.     b : vmatrixptr;
  819.     f : text;
  820.     titlestr: String;
  821. Begin
  822.     If Not FileExists( fid ) Then Begin 
  823.         writeln( 'trying to open : ', fid );
  824.         dispatch^.nrerror( 'reada: file does not exist' );
  825.     End;
  826.     
  827.     assign( f, fid );
  828.     reset( f );
  829.     
  830.     writeln( 'reading: ', fid );
  831.     readln( f, titlestr );
  832.     writeln( 'title:' );
  833.     writeln( titlestr );
  834.     read( f, i, j );
  835.     new( b, makematrix( i, j ) );
  836.     For i := 1 To b^.r Do
  837.         For j := 1 To b^.c Do
  838.             read( f, b^.mm( i, j )^ );
  839.     close( f );
  840.     dispatch^.push( b );
  841.     reada := dispatch^.returnmat;
  842. End; { reada }
  843.  
  844. Procedure writea;
  845. Var i,j,w,d  : integer;
  846.     f :  text;
  847. Begin
  848.     assign( f, fid );
  849.     rewrite( f );
  850.     
  851.     w := getwid;
  852.     d := getdec;
  853.     writeln( f, titlestr );
  854.     writeln( f, a^.r, ' ', a^.c );
  855.     For i := 1 To a^.r Do Begin 
  856.         For j := 1 To a^.c Do
  857.             write( f, a^.m( i, j ): w: d, ' ' );
  858.         writeln( f );
  859.     End;
  860.     close( f );
  861. End; { reada }
  862.  
  863. Var
  864.    printdec, printwid : integer;
  865.  
  866.  
  867.  
  868. Function getwid;
  869. Begin
  870.     getwid := printwid;
  871. End;
  872.  
  873. Function getdec;
  874. Begin
  875.     getdec := printdec;
  876. End;
  877.  
  878. Procedure setwid;
  879. Begin
  880.     If wid > 0 Then printwid := wid;
  881. End;
  882.  
  883. Procedure  setdec;
  884. Begin
  885.     If ( decimals > 0 ) And( decimals <= printwid ) Then
  886.         printdec := decimals;
  887. End;
  888.  
  889. Begin
  890.     setwid( 6 );
  891.     setdec( 2 );
  892.     new( dispatch, InitMatrixStack );
  893. End.
  894.  
  895.